library(ggplot2)
library(reshape2)
library(dplyr)
library(tidyr)
library(GGally)
library(grid)
"%&%" = function(a,b) paste(a,b,sep="")
source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan-Paper/scripts/Heather/make-figures/multiplot.R')
my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
out.dir <- '~/GitHub/GenArch/GenArchPaper/OTD_enrichment/'
Pull h2 estimate CI>0 genes for each tissue & rm .\d+ from ensid
for(tistype in c("Tissue-Specific","Tissue-Wide")){
tsfile <- my.dir %&% 'GTEx_' %&% tistype %&% '_local_h2_se_geneinfo_no_description.txt'
ts <- read.table(tsfile,sep='\t',header=T)
tislist <- scan(my.dir %&% 'GTEx_nine_tissues','c')
tislist <- gsub("-","",tislist)
tislist <- c("CrossTissue",tislist)
h2tislist <- "h2." %&% tislist
setislist <- "se." %&% tislist
geneinfo <- ts[,1:3]
for(i in seq_along(tislist)){
h2tis <- h2tislist[i]
setis <- setislist[i]
data <- ts %>% select(AssociatedGeneName, EnsemblGeneID, matches(h2tis), matches(setis)) %>% mutate(ensid=substr(EnsemblGeneID,1,15))
cidata <- data[data[,3]-data[,4]*2>0,] ##pull genes with non-zero confidence intervals
print(tislist[i] %&% " " %&% tistype %&% ": " %&% dim(cidata)[1] %&% " non-zero h2 CI genes")
write.table(cidata, file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_info.txt",quote=FALSE,row.names = FALSE)
write.table(cidata[,5], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_ensid_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
write.table(cidata[,1], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_non-zeroCIgenes_gene_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
for(thresh in c(0,0.01,0.05,0.1)){
h2data <- data[data[,3]>thresh,] ##pull genes with h2 > thresh
print(tislist[i] %&% " " %&% tistype %&% ": " %&% dim(h2data)[1] %&% " h2 > " %&% thresh %&% " genes")
write.table(h2data, file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_info.txt",quote=FALSE,row.names = FALSE)
write.table(h2data[,5], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_ensid_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
write.table(h2data[,1], file=out.dir %&% tislist[i] %&% "_" %&% tistype %&% "_h2_" %&% thresh %&% "genes_gene_list.txt",quote=FALSE,row.names = FALSE, col.names = FALSE)
}
}
}
## [1] "CrossTissue Tissue-Specific: 2938 non-zero h2 CI genes"
## [1] "CrossTissue Tissue-Specific: 17007 h2 > 0 genes"
## [1] "CrossTissue Tissue-Specific: 10282 h2 > 0.01 genes"
## [1] "CrossTissue Tissue-Specific: 5042 h2 > 0.05 genes"
## [1] "CrossTissue Tissue-Specific: 2715 h2 > 0.1 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 207 non-zero h2 CI genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 17007 h2 > 0 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 6908 h2 > 0.01 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 2463 h2 > 0.05 genes"
## [1] "AdiposeSubcutaneous Tissue-Specific: 650 h2 > 0.1 genes"
## [1] "ArteryTibial Tissue-Specific: 306 non-zero h2 CI genes"
## [1] "ArteryTibial Tissue-Specific: 17007 h2 > 0 genes"
## [1] "ArteryTibial Tissue-Specific: 7779 h2 > 0.01 genes"
## [1] "ArteryTibial Tissue-Specific: 3126 h2 > 0.05 genes"
## [1] "ArteryTibial Tissue-Specific: 892 h2 > 0.1 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 298 non-zero h2 CI genes"
## [1] "HeartLeftVentricle Tissue-Specific: 17007 h2 > 0 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 7796 h2 > 0.01 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 4012 h2 > 0.05 genes"
## [1] "HeartLeftVentricle Tissue-Specific: 1614 h2 > 0.1 genes"
## [1] "Lung Tissue-Specific: 223 non-zero h2 CI genes"
## [1] "Lung Tissue-Specific: 17007 h2 > 0 genes"
## [1] "Lung Tissue-Specific: 7143 h2 > 0.01 genes"
## [1] "Lung Tissue-Specific: 2738 h2 > 0.05 genes"
## [1] "Lung Tissue-Specific: 797 h2 > 0.1 genes"
## [1] "MuscleSkeletal Tissue-Specific: 434 non-zero h2 CI genes"
## [1] "MuscleSkeletal Tissue-Specific: 17007 h2 > 0 genes"
## [1] "MuscleSkeletal Tissue-Specific: 7371 h2 > 0.01 genes"
## [1] "MuscleSkeletal Tissue-Specific: 2356 h2 > 0.05 genes"
## [1] "MuscleSkeletal Tissue-Specific: 622 h2 > 0.1 genes"
## [1] "NerveTibial Tissue-Specific: 379 non-zero h2 CI genes"
## [1] "NerveTibial Tissue-Specific: 17006 h2 > 0 genes"
## [1] "NerveTibial Tissue-Specific: 7797 h2 > 0.01 genes"
## [1] "NerveTibial Tissue-Specific: 3491 h2 > 0.05 genes"
## [1] "NerveTibial Tissue-Specific: 1216 h2 > 0.1 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 449 non-zero h2 CI genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 17007 h2 > 0 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 8604 h2 > 0.01 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 3589 h2 > 0.05 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Specific: 1097 h2 > 0.1 genes"
## [1] "Thyroid Tissue-Specific: 442 non-zero h2 CI genes"
## [1] "Thyroid Tissue-Specific: 17007 h2 > 0 genes"
## [1] "Thyroid Tissue-Specific: 7757 h2 > 0.01 genes"
## [1] "Thyroid Tissue-Specific: 3141 h2 > 0.05 genes"
## [1] "Thyroid Tissue-Specific: 1010 h2 > 0.1 genes"
## [1] "WholeBlood Tissue-Specific: 490 non-zero h2 CI genes"
## [1] "WholeBlood Tissue-Specific: 17007 h2 > 0 genes"
## [1] "WholeBlood Tissue-Specific: 8504 h2 > 0.01 genes"
## [1] "WholeBlood Tissue-Specific: 2918 h2 > 0.05 genes"
## [1] "WholeBlood Tissue-Specific: 819 h2 > 0.1 genes"
## [1] "CrossTissue Tissue-Wide: 2938 non-zero h2 CI genes"
## [1] "CrossTissue Tissue-Wide: 17007 h2 > 0 genes"
## [1] "CrossTissue Tissue-Wide: 10282 h2 > 0.01 genes"
## [1] "CrossTissue Tissue-Wide: 5042 h2 > 0.05 genes"
## [1] "CrossTissue Tissue-Wide: 2715 h2 > 0.1 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 1127 non-zero h2 CI genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 17007 h2 > 0 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 6236 h2 > 0.01 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 3407 h2 > 0.05 genes"
## [1] "AdiposeSubcutaneous Tissue-Wide: 1736 h2 > 0.1 genes"
## [1] "ArteryTibial Tissue-Wide: 1171 non-zero h2 CI genes"
## [1] "ArteryTibial Tissue-Wide: 17007 h2 > 0 genes"
## [1] "ArteryTibial Tissue-Wide: 6314 h2 > 0.01 genes"
## [1] "ArteryTibial Tissue-Wide: 3546 h2 > 0.05 genes"
## [1] "ArteryTibial Tissue-Wide: 1889 h2 > 0.1 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 749 non-zero h2 CI genes"
## [1] "HeartLeftVentricle Tissue-Wide: 17007 h2 > 0 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 6352 h2 > 0.01 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 3865 h2 > 0.05 genes"
## [1] "HeartLeftVentricle Tissue-Wide: 2179 h2 > 0.1 genes"
## [1] "Lung Tissue-Wide: 929 non-zero h2 CI genes"
## [1] "Lung Tissue-Wide: 17007 h2 > 0 genes"
## [1] "Lung Tissue-Wide: 6331 h2 > 0.01 genes"
## [1] "Lung Tissue-Wide: 3387 h2 > 0.05 genes"
## [1] "Lung Tissue-Wide: 1693 h2 > 0.1 genes"
## [1] "MuscleSkeletal Tissue-Wide: 1063 non-zero h2 CI genes"
## [1] "MuscleSkeletal Tissue-Wide: 17007 h2 > 0 genes"
## [1] "MuscleSkeletal Tissue-Wide: 6049 h2 > 0.01 genes"
## [1] "MuscleSkeletal Tissue-Wide: 2814 h2 > 0.05 genes"
## [1] "MuscleSkeletal Tissue-Wide: 1331 h2 > 0.1 genes"
## [1] "NerveTibial Tissue-Wide: 1466 non-zero h2 CI genes"
## [1] "NerveTibial Tissue-Wide: 17006 h2 > 0 genes"
## [1] "NerveTibial Tissue-Wide: 6594 h2 > 0.01 genes"
## [1] "NerveTibial Tissue-Wide: 4057 h2 > 0.05 genes"
## [1] "NerveTibial Tissue-Wide: 2399 h2 > 0.1 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 1198 non-zero h2 CI genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 17007 h2 > 0 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 6504 h2 > 0.01 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 3530 h2 > 0.05 genes"
## [1] "SkinSunExposed(Lowerleg) Tissue-Wide: 1809 h2 > 0.1 genes"
## [1] "Thyroid Tissue-Wide: 1327 non-zero h2 CI genes"
## [1] "Thyroid Tissue-Wide: 17006 h2 > 0 genes"
## [1] "Thyroid Tissue-Wide: 6590 h2 > 0.01 genes"
## [1] "Thyroid Tissue-Wide: 3758 h2 > 0.05 genes"
## [1] "Thyroid Tissue-Wide: 2082 h2 > 0.1 genes"
## [1] "WholeBlood Tissue-Wide: 945 non-zero h2 CI genes"
## [1] "WholeBlood Tissue-Wide: 17007 h2 > 0 genes"
## [1] "WholeBlood Tissue-Wide: 5426 h2 > 0.01 genes"
## [1] "WholeBlood Tissue-Wide: 2586 h2 > 0.05 genes"
## [1] "WholeBlood Tissue-Wide: 1261 h2 > 0.1 genes"
write(data$ensid, file=out.dir %&% "full_tested_ensid_list.txt",ncolumns=1)
Define known genes for 7 WTCCC diseases based on the GWAS catalog and make list of ALL genes in catalog
#catalog used:
gwasfile <- my.dir %&% 'gwas_catalog_v1.0-downloaded_2015-06-02.tsv'
# gwas <- data.frame(read.table(gwasfile,header=TRUE,sep='\t',quote="",comment.char="",as.is=TRUE)) #for reference
runPERL <- "perl " %&% my.dir %&% "24_define_disease_genes.pl " %&% gwasfile
system(runPERL)
Disease gene enrichment in top h2 genes by tissue
lci<-function(x) quantile(x, c(.025, 0.975),na.rm=T)[[1]]
uci<-function(x) quantile(x, c(.025, 0.975),na.rm=T)[[2]]
##is testvec signficanty enriched for setvec? make n samples of size(testvec) from fullvec and count overlap
enrichment <- function(setvec, testvec, fullvec, nperms = 1000){
obs <- length(testvec[testvec %in% setvec])
counts <- vector()
for(i in 1:nperms){
rantest <- base::sample(fullvec,length(testvec),replace=FALSE)
cnt <- length(rantest[rantest %in% setvec])
counts <- c(counts,cnt)
}
empp <- mean(counts>obs)
meanc <- mean(counts)
lc <- lci(counts)
uc <- uci(counts)
return(c(obs,meanc,lc,uc,empp))
}
set.seed(42)
fullgenelist <- as.character(ts$AssociatedGeneName)
dislist <- c("BD","CAD","HT","T1D","T2D","CD","RA","ALL")
nperms <- 10000
for(thresh in c("non-zeroCI","h2_0.1","h2_0.05","h2_0.01")){
results <- data.frame(type=character(0),dis=character(0),tis=character(0),obsOverlap=double(0),meanOverlap=double(0),lCI=double(0),uCI=double(0),empPval=double(0))
for(tistype in c("Tissue-Specific","Tissue-Wide")){
for(i in seq_along(dislist)){
dis <- dislist[i]
disgenes <- scan(out.dir %&% "gwas." %&% dis %&% ".tsv","character")
for(j in seq_along(tislist)){
tis <- tislist[j]
tisgenes <- scan(out.dir %&% tis %&% '_' %&% tistype %&% '_' %&% thresh %&% 'genes_gene_list.txt','character')
res <- enrichment(disgenes,tisgenes,fullgenelist,nperms=nperms)
resvec <- data.frame(type=tistype,dis=dis,tis=tis,obsOverlap=res[1],meanOverlap=res[2],lCI=res[3],uCI=res[4],empPval=res[5])
results <- rbind(results,resvec)
}
}
}
sortres <- arrange(results,empPval)
write.table(sortres,file=out.dir %&% "GWAS_catalog_disease_gene_enrichment_in_" %&% thresh %&% "_genes_by_tissue_" %&% nperms %&% "perms.txt",row.names=FALSE,quote=FALSE)
print("GWAS catalog disease gene enrichment in " %&% thresh %&% " genes by tissue " %&% nperms %&% " perms:")
print(head(sortres,n=20L))
}
## [1] "GWAS catalog disease gene enrichment in non-zeroCI genes by tissue 10000 perms:"
## type dis tis obsOverlap meanOverlap lCI uCI
## 1 Tissue-Specific ALL Lung 129 102.1142 88 117
## 2 Tissue-Specific T2D Lung 9 2.1773 0 6
## 3 Tissue-Specific CD MuscleSkeletal 19 9.5357 4 16
## 4 Tissue-Specific BD Thyroid 17 8.2913 3 14
## 5 Tissue-Specific ALL ArteryTibial 165 140.1439 123 157
## 6 Tissue-Specific HT Lung 3 0.5509 0 2
## 7 Tissue-Specific ALL HeartLeftVentricle 160 136.3636 120 153
## 8 Tissue-Specific T2D AdiposeSubcutaneous 6 2.0594 0 5
## 9 Tissue-Specific T1D WholeBlood 7 2.7214 0 6
## 10 Tissue-Specific HT HeartLeftVentricle 3 0.7396 0 3
## 11 Tissue-Specific T2D HeartLeftVentricle 7 2.9590 0 7
## 12 Tissue-Specific CAD HeartLeftVentricle 8 3.5693 1 8
## 13 Tissue-Wide BD Thyroid 36 24.9775 16 35
## 14 Tissue-Specific T2D Thyroid 8 4.3519 1 9
## 15 Tissue-Specific RA Lung 5 2.3356 0 6
## 16 Tissue-Specific RA HeartLeftVentricle 6 3.1070 0 7
## 17 Tissue-Specific ALL Thyroid 219 202.3259 182 223
## 18 Tissue-Specific BD ArteryTibial 9 5.7380 2 11
## 19 Tissue-Specific RA WholeBlood 8 5.1027 1 10
## 20 Tissue-Specific T1D HeartLeftVentricle 3 1.6345 0 4
## empPval
## 1 0.0001
## 2 0.0002
## 3 0.0014
## 4 0.0015
## 5 0.0015
## 6 0.0024
## 7 0.0027
## 8 0.0049
## 9 0.0054
## 10 0.0055
## 11 0.0111
## 12 0.0120
## 13 0.0120
## 14 0.0275
## 15 0.0331
## 16 0.0402
## 17 0.0478
## 18 0.0643
## 19 0.0708
## 20 0.0794
## [1] "GWAS catalog disease gene enrichment in h2_0.1 genes by tissue 10000 perms:"
## type dis tis obsOverlap meanOverlap
## 1 Tissue-Specific HT HeartLeftVentricle 11 3.9737
## 2 Tissue-Specific BD ArteryTibial 28 16.7268
## 3 Tissue-Specific BD Thyroid 30 18.9886
## 4 Tissue-Specific T2D Thyroid 18 9.9539
## 5 Tissue-Specific CD MuscleSkeletal 23 13.6412
## 6 Tissue-Specific ALL Thyroid 498 462.4593
## 7 Tissue-Specific ALL Lung 396 365.0230
## 8 Tissue-Specific T2D AdiposeSubcutaneous 12 6.3959
## 9 Tissue-Specific T1D WholeBlood 9 4.5149
## 10 Tissue-Specific T2D SkinSunExposed(Lowerleg) 17 10.8103
## 11 Tissue-Wide BD Thyroid 50 39.2137
## 12 Tissue-Specific ALL ArteryTibial 433 408.3995
## 13 Tissue-Specific HT Lung 4 1.9382
## 14 Tissue-Specific ALL MuscleSkeletal 304 284.5559
## 15 Tissue-Specific CAD AdiposeSubcutaneous 12 7.8481
## 16 Tissue-Specific ALL AdiposeSubcutaneous 317 297.4614
## 17 Tissue-Specific T2D NerveTibial 17 12.0049
## 18 Tissue-Specific BD AdiposeSubcutaneous 17 12.1849
## 19 Tissue-Wide HT Thyroid 8 5.1605
## 20 Tissue-Wide BD SkinSunExposed(Lowerleg) 42 34.1768
## lCI uCI empPval
## 1 1.000 8 0.0005
## 2 9.000 25 0.0036
## 3 11.000 27 0.0051
## 4 4.000 16 0.0052
## 5 7.000 21 0.0061
## 6 432.000 492 0.0114
## 7 338.000 392 0.0119
## 8 2.000 12 0.0123
## 9 1.000 9 0.0146
## 10 5.000 17 0.0248
## 11 28.000 51 0.0307
## 12 380.000 437 0.0401
## 13 0.000 5 0.0408
## 14 261.000 308 0.0473
## 15 3.000 14 0.0517
## 16 273.975 322 0.0544
## 17 6.000 19 0.0547
## 18 6.000 19 0.0643
## 19 1.000 10 0.0669
## 20 24.000 45 0.0685
## [1] "GWAS catalog disease gene enrichment in h2_0.05 genes by tissue 10000 perms:"
## type dis tis obsOverlap meanOverlap
## 1 Tissue-Specific HT Lung 15 6.7802
## 2 Tissue-Wide T2D ArteryTibial 52 34.9891
## 3 Tissue-Specific ALL ArteryTibial 1507 1430.8166
## 4 Tissue-Specific ALL Thyroid 1500 1438.3020
## 5 Tissue-Wide HT ArteryTibial 15 8.8099
## 6 Tissue-Specific BD NerveTibial 83 65.6931
## 7 Tissue-Specific ALL Lung 1305 1253.4985
## 8 Tissue-Specific CD ArteryTibial 83 68.4924
## 9 Tissue-Wide ALL MuscleSkeletal 1333 1287.7332
## 10 Tissue-Specific T2D Lung 35 27.0397
## 11 Tissue-Specific CAD AdiposeSubcutaneous 38 29.7497
## 12 Tissue-Wide ALL ArteryTibial 1669 1623.8195
## 13 Tissue-Specific HT HeartLeftVentricle 14 9.8989
## 14 Tissue-Specific ALL MuscleSkeletal 1115 1078.7831
## 15 Tissue-Specific ALL HeartLeftVentricle 1883 1836.9226
## 16 Tissue-Specific HT MuscleSkeletal 9 5.8051
## 17 Tissue-Wide HT Thyroid 13 9.3239
## 18 Tissue-Specific ALL SkinSunExposed(Lowerleg) 1682 1642.5537
## 19 Tissue-Specific T2D Thyroid 38 30.9843
## 20 Tissue-Specific ALL NerveTibial 1637 1598.3363
## lCI uCI empPval
## 1 3 12 0.0004
## 2 25 46 0.0008
## 3 1381 1480 0.0010
## 4 1389 1486 0.0051
## 5 4 14 0.0077
## 6 52 80 0.0082
## 7 1207 1300 0.0154
## 8 55 83 0.0238
## 9 1240 1336 0.0323
## 10 18 37 0.0395
## 11 20 40 0.0418
## 12 1572 1676 0.0434
## 13 5 15 0.0480
## 14 1034 1122 0.0499
## 15 1782 1892 0.0501
## 16 2 10 0.0601
## 17 4 15 0.0635
## 18 1591 1695 0.0666
## 19 22 41 0.0686
## 20 1547 1650 0.0687
## [1] "GWAS catalog disease gene enrichment in h2_0.01 genes by tissue 10000 perms:"
## type dis tis obsOverlap meanOverlap lCI
## 1 Tissue-Specific ALL ArteryTibial 3670 3560.6350 3498.000
## 2 Tissue-Specific HT HeartLeftVentricle 27 19.2386 13.000
## 3 Tissue-Wide ALL ArteryTibial 2961 2889.8261 2829.000
## 4 Tissue-Specific ALL NerveTibial 3642 3569.2810 3507.000
## 5 Tissue-Wide BD NerveTibial 141 124.0499 107.000
## 6 Tissue-Specific BD MuscleSkeletal 156 138.8364 122.000
## 7 Tissue-Wide T2D ArteryTibial 74 62.3432 50.000
## 8 Tissue-Wide HT ArteryTibial 21 15.5968 10.000
## 9 Tissue-Specific ALL HeartLeftVentricle 3626 3568.2826 3506.000
## 10 Tissue-Specific ALL Lung 3326 3270.1042 3205.975
## 11 Tissue-Specific HT NerveTibial 24 19.2605 13.000
## 12 Tissue-Wide HT Thyroid 21 16.3225 10.000
## 13 Tissue-Wide BD AdiposeSubcutaneous 130 117.2681 101.000
## 14 Tissue-Specific ALL Thyroid 3600 3550.8891 3487.000
## 15 Tissue-Wide T2D Thyroid 74 65.0499 53.000
## 16 Tissue-Specific ALL AdiposeSubcutaneous 3207 3161.6611 3099.000
## 17 Tissue-Specific ALL MuscleSkeletal 3419 3374.3486 3311.000
## 18 Tissue-Wide CAD MuscleSkeletal 82 72.9027 60.000
## 19 Tissue-Specific CD AdiposeSubcutaneous 164 151.5080 133.000
## 20 Tissue-Wide CD WholeBlood 130 118.8994 102.000
## uCI empPval
## 1 3625.000 0.0003
## 2 26.000 0.0052
## 3 2952.000 0.0108
## 4 3633.000 0.0140
## 5 141.000 0.0229
## 6 156.000 0.0233
## 7 75.000 0.0286
## 8 22.000 0.0301
## 9 3631.000 0.0350
## 10 3333.000 0.0379
## 11 26.000 0.0501
## 12 23.000 0.0509
## 13 134.000 0.0614
## 14 3615.025 0.0646
## 15 78.000 0.0692
## 16 3224.000 0.0740
## 17 3438.000 0.0812
## 18 86.000 0.0822
## 19 170.000 0.0852
## 20 136.000 0.0960
distribution of h2 for disease vs non disease
dislist <- c("BD","CAD","HT","T1D","T2D","CD","RA","ALL")
tislist <- c("CrossTissue","AdiposeSubcutaneous","ArteryTibial","HeartLeftVentricle","Lung","MuscleSkeletal","NerveTibial","SkinSunExposed(Lowerleg)","Thyroid","WholeBlood")
typelist<-c("Tissue-Specific","Tissue-Wide")
for(thresh in c(0.1,0.05,0)){
for(tistype in typelist){
for(tis in tislist){
info <- read.table(out.dir %&% tis %&% '_' %&% tistype %&% '_h2_' %&% thresh %&% 'genes_info.txt',header=T)
finaldf <- data.frame(AssociatedGeneName=character(0),EnsemblGeneID=character(0),h2=double(0),se=double(0),ensid=character(0),diseaseGene=logical(0L),disease=character(0))
for(dis in dislist){
setvec <- scan(out.dir %&% "gwas." %&% dis %&% ".tsv","character")
disinfo <- info %>% mutate(diseaseGene=(info[,1] %in% setvec),disease=dis)
colnames(disinfo) <- c("AssociatedGeneName","EnsemblGeneID","h2","se","ensid","diseaseGene","disease")
finaldf <- rbind(finaldf,disinfo)
}
p<-ggplot(finaldf,aes(x=finaldf[,3],fill=diseaseGene,color=diseaseGene)) + facet_wrap(~disease,ncol=2) + geom_density(alpha=0.3) + xlab("h2") + ggtitle(tistype %&% ' ' %&% tis %&% ' h2 > ' %&% thresh)
print(p)
p<-ggplot(finaldf,aes(y=finaldf[,3],x=diseaseGene)) + facet_wrap(~disease,ncol=4) + geom_jitter(aes(colour=diseaseGene),alpha=0.3,position = position_jitter(width = .15)) + geom_boxplot(alpha=0,outlier.size=NA) + xlab("diseaseGene") + ylab("h2") + ggtitle(tistype %&% ' ' %&% tis %&% ' h2 > ' %&% thresh) +theme_bw() + theme(legend.position="none")
print(p)
}
}
}
